
function Probm = update_dist_int_edu_mto(myDist, Pr, Ss, rs, tax, w_mto, Hs, flag)
    global bet0 bet1 tau a b w e pi_theta N_a N_theta N_w N_e w_max w_min alpha xi xi1 ub lb pi_e...
        w theta Int_pi_w sigmae elasts log_e_mu log_e_sigma phi sigma gamma rts w_mult wspread unemp skill_shock
    
    % Pre-allocate pdf matrix and caluclate variables
    Probm = zeros(size(myDist));
    
    pi_w = sum(myDist, 2:3);
    w_bar = dot(w, pi_w ); 

    b0 = b;

    % Choice Probabilities
    [~, tax, mean_wage] = find_rents_mto(rs, Ss,  Hs, myDist, w_mto);
    [P_S, F_es, wprimeS, ~] = choiceProb_ct_par(Ss, rs, tax, w_mto, mean_wage);
    
    % Update parents' distribution
    if flag == 1
        F_es = F_es*skill_shock;
    end

   
    myDist1  = sum(bsxfun(@times, myDist .* P_S , reshape(pi_theta, [1,1,1,1,2])),5);
    w_bar_p = sum(myDist1 .* wprimeS, 1:4);
 
    w2 = w_bar_p * wspread;
    w_max = max(w2);
    w_min = min(w2);
    ln_e_hat = -1/2*log_e_sigma^2;
    
    % Update distribution
    for c1 = 1:length(Ss)
        for t = 1:2
            for j = 1:N_a
                for i = 1:N_w
                    wage = w(i);               
                    t1 = unemp + (b   + squeeze(F_es(i,j,:,c1))') .* f_w(wage*(1-tax));
                    midpts = log(w2) - log(t1);
                    pts = ([midpts(1,:);(midpts(2:end,:) + midpts(1:end-1,:))/2; Inf Inf Inf] - ln_e_hat)/log_e_sigma;
                    mean_prob = diff(normcdf(pts)) .* squeeze(myDist(i, j, :))';
                    tmp = bsxfun(@times, mean_prob*squeeze(P_S(i,j,:,c1,t)), Pr(j,:));
                    pr1 =  tmp .* pi_theta(t);
                    Probm(:, :, c1) = Probm(:, :, c1) + pr1;

                end
            end
        end
    end
    w = w2;
end